% For a DSGE model approximated up to third order, this function calculates 
% the unconditional first and second moments based on the the pruning
% scheme. The system reads
%     xf(t+1)  = hx*xf(t)  + sig*eta*eps(t);
%     xs(t+1)  = hx*xs(t)  + 0.5*reshape(hxx,nx,nx^2)*kron(xf(t),xf(t))+ 1/2*sig^2*hss;
%     xrd(t+1) = hx*xrd(t) + 2*0.5*reshape(hxx,nx,nx^2)*kron(xf(t),xs(t))...
%                + 1/6*reshape(hxxx,nx,nx^3)*kron(xf(t),kron(xf(t),xf(t)))
%                + 3/6*hssx*sig^2*xf(t) + 1/6*sig^3*hsss;
%     y(t+1)  = gx*(xf(t+1)+xs(t+1)+xrd(t+1)) ...
%                + 0.5*reshape(gxx,ny,nx^2)*(kron(xf(t+1),xf(t+1)+2*kron(xf(t+1),xs(t+1))+ 1/2*sig^2*gss
%                + 1/6*reshape(gxxx,ny,nx^3)*(kron(xf(t+1),kron(xf(t+1),xf(t+1)))
%                + 3/6*gssx*sig^2*xf(t+1) + 1/6*sig^3*gsss;
%
% IMPORTANT: Here we allow for general shock distributions, i.e. they may be
% non-symmetric. Here the Lyaponov method is used to solve for the
% variances
%
% INPUTS:
% gx,...,sig : The second order approximation of the DSGE model
% vectorMom3 : vector with dimensions 1*ne with third moments of shocks
% vectorMom4 : vector with dimensions 1*ne with fourth moments of shocks
% vectorMom5 : vector with dimensions 1*ne with fifth moments of shocks
% vectorMom6 : vector with dimensions 1*ne with sixth moments of shocks
% auto_lag   : the lag lenght for the auto-covariances (they need to be divided
%              by the standard deviations to get the correlations)
%
% OUTPUTS: 
% See below
%
% By M. Andreasen, 2015

function out = UnconditionalMoments_3rd_Lyap(gx,gxx,gss,gxxx,gssx,gsss,hx,hxx,hss,hxxx,hssx,hsss,eta,sig,...
      vectorMom3,vectorMom4,vectorMom5,vectorMom6,Var_xs,Var_xfxf,Var_xs_xfxf,auto_lag,Var_zold);

% The dimensions
ny      = size(gx,1);
[nx,ne] = size(eta);
sigeta  = sig*eta;

% The value of E_u
E_eps3 = zeros(ne^3,1);
index = 0;
for phi1=1:ne
    for phi2=1:ne
        for phi3=1:ne
            index = index + 1;
            if phi1 == phi2 && phi1 == phi3
                E_eps3(index,1) = vectorMom3(1,phi1);
            end
        end
    end
end
E_u = kron(sigeta,kron(sigeta,sigeta))*E_eps3;

% The auxiliary system z=[xf xs kron(xf,xf) xrd kron(xf,xs) kron(xf,kron(xf,xf))]
% The loading matrix
A = [hx                    ,zeros(nx,nx)    ,zeros(nx,nx^2)          , zeros(nx,nx)  ,  zeros(nx,nx^2)            , zeros(nx,nx^3)  ;
    zeros(nx,nx)          ,hx              ,0.5*reshape(hxx,nx,nx^2), zeros(nx,nx)  ,  zeros(nx,nx^2)            , zeros(nx,nx^3)  ;
    zeros(nx^2,nx)        ,zeros(nx^2,nx)  ,kron(hx,hx)             , zeros(nx^2,nx),  zeros(nx^2,nx^2)          , zeros(nx^2,nx^3);
    3/6*hssx*sig^2        ,zeros(nx,nx)    ,zeros(nx,nx^2)          , hx            ,  2*0.5*reshape(hxx,nx,nx^2), 1/6*reshape(hxxx,nx,nx^3);
    kron(hx,0.5*hss*sig^2),zeros(nx^2,nx)  ,zeros(nx^2,nx^2)        , zeros(nx^2,nx),  kron(hx,hx)               , kron(hx,0.5*reshape(hxx,nx,nx^2));
    zeros(nx^3,nx)        ,zeros(nx^3,nx)  ,zeros(nx^3,nx^2)        , zeros(nx^3,nx),  zeros(nx^3,nx^2)          , kron(hx,kron(hx,hx))   ];
% The constant term
c = [zeros(nx,1);
    0.5*hss*sig^2;
    kron(sig*eta,sig*eta)*reshape(eye(ne),ne^2,1);
    1/6*hsss*sig^3;
    zeros(nx^2,1);
    E_u];

% the mean value in the auxiliary system
%E_z       = inv(sparse(eye(3*nx+2*nx^2+nx^3)-A))*c;
E_z        = (sparse((eye(3*nx+2*nx^2+nx^3)-A)))\c;
E_xs       = E_z(nx+1:2*nx,1);
E_xf_xf    = reshape(E_z(2*nx+1:2*nx+nx^2,1),nx,nx);
E_xrd      = E_z(2*nx+nx^2+1:3*nx+nx^2,1);
E_xf_xs    = reshape(E_z(3*nx+nx^2+1:3*nx+2*nx^2,1),nx,nx);
E_xf_xf_xf = reshape(E_z(3*nx+2*nx^2+1:3*nx+2*nx^2+nx^3,1),nx,nx,nx);

% Computing variance of the innovations
Var_inov = Get_VarInov_3rd(hx,hxx,hss,sigeta,sig,...
    vectorMom3,vectorMom4,vectorMom5,vectorMom6,E_xs,E_xf_xf,E_xf_xs,E_xf_xf_xf,E_u,...
    Var_xs,Var_xfxf,Var_xs_xfxf,nx,ne);


% Computing B*cov(z_t,xi_t+1+s)
LeadS = 0;
%BCov_xiLeadS_z = Get_BCov_xiLeadS_z(hx,hxx,hxxx,hssx,sigeta,sig,...
%    E_xs,E_xf_xf,E_xf_xs,E_xf_xf_xf,Var_xfxf,Var_xs_xfxf,nx,ne,LeadS);
BCov_xiLeadS_z = Get_BCov_xiLeadS_z_v2(hx,hxx,hxxx,hssx,sigeta,sig,...
    E_xs,E_xf_xf,E_xf_xs,E_xf_xf_xf,Var_xfxf,Var_xs_xfxf,nx,ne,LeadS);

% Constant terms
const_term = Var_inov + BCov_xiLeadS_z*A' + A*BCov_xiLeadS_z';

% Computing the variances of the state variables, using the Lyaponov routine
Var_z      = dlyap_doubling(A,const_term,Var_zold);
Var_xf     = Var_z(1:nx,1:nx);
Var_xs     = Var_z(nx+1:2*nx,nx+1:2*nx);
Var_xrd    = Var_z(2*nx+nx^2+1:3*nx+nx^2,2*nx+nx^2+1:3*nx+nx^2);
Var_x      = Var_z(1:nx,1:nx) + Var_z(nx+1:2*nx,nx+1:2*nx) + Var_xrd ...
             + Var_z(1:nx,nx+1:2*nx)...                   %Cov(xf,xs)
             + Var_z(nx+1:2*nx,1:nx)...                   %Cov(xs,xf)
             + Var_z(1:nx,2*nx+nx^2+1:3*nx+nx^2)...       %Cov(xf,xrd)
             + Var_z(2*nx+nx^2+1:3*nx+nx^2,1:nx)...       %Cov(xrd,xf)             
             + Var_z(nx+1:2*nx,2*nx+nx^2+1:3*nx+nx^2)...  %Cov(xs,xrd)
             + Var_z(2*nx+nx^2+1:3*nx+nx^2,nx+1:2*nx);    %Cov(xrd,xs)
        
% Computing the auto-covariances
Cov_z      = zeros(size(Var_z,1),size(Var_z,2),auto_lag);
Cov_x      = zeros(nx,nx,auto_lag);
for i=1:auto_lag
    if i == 1
        Cov_z(:,:,i) = A*Var_z+BCov_xiLeadS_z;
    else
        LeadS = i-1;
        BCov_xiLeadS_z = Get_BCov_xiLeadS_z_v2(hx,hxx,hxxx,hssx,sigeta,sig,...
            E_xs,E_xf_xf,E_xf_xs,E_xf_xf_xf,Var_xfxf,Var_xs_xfxf,nx,ne,LeadS);
        Cov_z(:,:,i) = A*Cov_z(:,:,i-1)+BCov_xiLeadS_z;
    end
    Cov_xf      = Cov_z(1:nx,1:nx,i);
    Cov_xs      = Cov_z(nx+1:2*nx,nx+1:2*nx,i);
    Cov_xrd     = Cov_z(2*nx+nx^2+1:3*nx+nx^2,2*nx+nx^2+1:3*nx+nx^2,i);
    Cov_x(:,:,i)= Cov_xf + Cov_xs + Cov_xrd ...
        + Cov_z(1:nx,nx+1:2*nx)...                   %Cov(xf,xs)
        + Cov_z(nx+1:2*nx,1:nx)...                   %Cov(xs,xf)
        + Cov_z(1:nx,2*nx+nx^2+1:3*nx+nx^2)...       %Cov(xf,xrd)
        + Cov_z(2*nx+nx^2+1:3*nx+nx^2,1:nx)...       %Cov(xrd,xf)
        + Cov_z(nx+1:2*nx,2*nx+nx^2+1:3*nx+nx^2)...  %Cov(xs,xrd)
        + Cov_z(2*nx+nx^2+1:3*nx+nx^2,nx+1:2*nx);    %Cov(xrd,xs)   
end

% Moments of y: Recall y = D*z+d
D          = [gx+3/6*gssx*sig^2 gx 0.5*reshape(gxx,ny,nx^2) gx 2*0.5*reshape(gxx,ny,nx^2) 1/6*reshape(gxxx,ny,nx^3)];
d          = 0.5*gss*sig^2+1/6*gsss*sig^3;
E_y        = D*E_z+d;
Var_y      = D*Var_z*D';
Cov_y      = zeros(ny,ny,auto_lag);
for i=1:auto_lag
    Cov_y(:,:,i) = D*Cov_z(:,:,i)*D';
end

% output
out.E_z     = E_z;
out.E_xrd   = E_xrd;
out.Var_z   = Var_z;
out.Var_x   = Var_x;
out.Cov_z   = Cov_z;
out.Cov_x   = Cov_x;
out.E_y     = E_y;
out.Var_y   = Var_y;
out.Cov_y   = Cov_y;

end



